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A recently introduced family of lattice Boltzmann (LB) models (Karlin, Bosch, Chikatamarla, 
Phys. Rev. E, 2014; Ref|22)l is studied in detail for incompressible two-dimensional flows. A 
framework for developing LB models based on entropy considerations is laid out extensively. Second 
order rate of convergence is numerically confirmed and it is demonstrated that these entropy based 
models recover the Navier-Stokes solution in the hydrodynamic limit. Comparison with the standard 
Bhatnagar-Gross-Krook (LBGK) and the entropic lattice Boltzmann method (ELBM) demonstrates 
the superior stability and accuracy for several benchmark flows and a range of grid resolutions 
and Reynolds numbers. High Reynolds number regimes are investigated through the simulation 
of two-dimensional turbulence, particularly for under-resolved cases. Compared to resolved LBGK 
simulations, the presented class of LB models demonstrate excellent performance and capture the 
turbulence statistics with good accuracy. 


I. INTRODUCTION 

In recent years, the lattice Boltzmann (LB) method has made substantial progress towards a 
successful and particularly efficient approach to computational fluid dynamics. By employing a 
system of discrete kinetic equations rather than solving directly the macroscopic flow equations it 
has shown its potential in a wide range of applications, from turbulence phenomena [SI [U to 
flows at a micron scale [5] and multiphase flows [33j . to relativistic hydrodynamics m, soft-glassy 
systems |5] and beyond. 

The kinetic system describes the discrete-time dynamics of populations which are 

designed to reproduce the Navier-Stokes equations in the hydrodynamic limit. Each population fi 
corresponds to a discrete microscopic velocity Vi, i = 1,... ,b, which fits into a regular spatial lattice 
with the nodes x. This enables a highly efficient ‘stream-along-links-and-equilibrate-at-nodes’ 
realization of the LB algorithm. We consider the single-phase isothermal case in two dimensions 
for the purpose or this paper. A general form of the LB equation can be written as 

f,{x + Vi,t+l) = f- = (l-/3)/,(a;,t) +/3/™(a;,t). (1) 

Here the left-hand side is the propagation of the populations along the lattice links, while the 
right-hand side is the so-called post-collision state /'. The post-collision state is a convex linear 
combination between the pre-collision state / and the mirror state The choice of as 

/r’' = 2/r-/z (2) 

results in the well known LBGK model which has most notably led to the success of the method, 
in particular for the simulation of incompressible flows ladniiniEn] and complex hydrodynamic 
phenomena |T115^. 

The local equilibrium is found as a maximizer of the entropy, 

sm-gP) 
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subject to fixed locally conserved fields, p = Y}i=i fi (density) and pu = 'Z,i=i Vifi (momentum 
density), and where the weights Wi are lattice-specific constants. The equilibrium can be 
approximated in closed form to second order in velocity as 


fr=w.p[ 


1 + 




‘i-cis. 


2ci 


+0(u^) 


( 4 ) 


where Cg is the speed of sound (a lattice dependent 0(1) constant). The LBGK equations, Q 
and (§> recover the Navier-Stokes equation for the fluid velocity u in the hydrodynamic limit, 
under the assumption u « Cg, while the relaxation parameter /3 € [0,1] is defined by the kinematic 
viscosity v 


V = 




( 5 ) 


In order to achieve high Reynolds numbers, one must decrease viscosity as the flow velocity 
is restricted to low Mach numbers by construction of the kinetic system and, thus, the limit 
/3 ^ 1 is of great importance. While the time step and the lattice spacing are connected by the 
relation 8xi = btvi (where 6xi is the lattice spacing in the direction of discrete velocity Vi), the 
kinematic viscosity can be chosen independently of the spatial and temporal resolution. This 
feature makes the LBGK model particularly attractive for high Reynolds number simulations, if 
only in principle. 

Despite of its promising nature and popularity, however, the LBGK model shows numerical 
instabilities already at moderate Reynolds numbers unless a rather high resolution is employed, 
which quickly becomes computationally prohibitive. This precluded the LB method from making 
a sustainable impact in the field of computational fluid dynamics. 

A number of approaches can be found in the literature intended to alleviate this issue. We 
will restrict the following short discussion to methods without explicit turbulence models. Most 
notably, the entropic lattice Boltzmann method (ELBM) features non-linear stability and has 
shown excellent performance [Hiisiisg. While ELBM converges to LBGK in the resolved case, 
it locally modifies the relaxation rate which in turn modifies the viscosity in order to fulfil the 
second law of thermodynamics by both enhancing and smoothing the features of the flow where 
necessary subject to an entropy condition. Another approach using multiple relaxation parameters 
(MRT) is widely used, which does not affect viscosity in the first place but requires careful tuning 
of the relaxation parameters. Although MRT models were successful in slightly stabilizing the 
LB method, they still remain challenged by high Reynolds numbers m- 

Recently, the authors have developed a scheme without the need for tunable parameters 
or turbulent viscosity (Karlin, Bosch, Ghikatamarla, Phys. Rev. E 2014; Ref |22] 1 which has 
demonstrated a significant extension in the operation range for simulations at high Reynolds 
numbers. Promising results have been reported for both two and three dimensions, as well as for 
complex boundaries and in presence of turbulence. Much alike ELBM, entropic considerations 
have been employed to render the scheme stable without introducing considerable computational 
overhead and by keeping the simplicity and locality of the LBGK and MRT schemes. Below we 
shall refer to this class of models as KBG models for brevity. 

While in |22j one particular realization of the model was discussed, this paper focuses on 
four variations of the KBG family which will be investigated in detail, both numerically and 
analytically. We restrict ourselves to two dimensional fully periodic domains in the absence of 
wall boundaries in order to assess the stability and accuracy of the scheme independently of the 
errors arising from the wall. 

The MRT class of LB models separate the relaxation into various groups based on separation 
of scales between the fast and slow varying moments. Moreover, since the relaxation of the 
off-diagonal parts of the pressure tensor are fixed by the choice of kinematic viscosity, the MRT 
scheme asserts that the relaxation of higher order moments should not affect the flow field (up to 
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the Navier-Stokes level) and hence can be used to construct more stable LB schemes. Following 
this line of thinking, several MRT schemes were suggested for the choice of relaxation of higher 
order moments (beyond the pressure tensor) [12] UHl [20] ■ The KBC models extend this idea 
using local entropy considerations and demonstrate that much higher Reynolds numbers can be 
achieved on much smaller grid sizes. 

Three well-studied benchmark flows - Green-Taylor vortex, doubly periodic shear layer and 
decaying two-dimensional turbulence - are simulated for all the four KBC models as well as for 
the classical LBGK and the non-linearly stable ELBM. Results are compared to each other as 
well as to reference solutions in order to evaluate the models. 


II. MOMENT REPRESENTATION 

We consider the standard nine-velocity model (D2Q9). The discrete velocities are constructed 
as a tensor product of two one-dimensional velocity sets, = i, where i = 0,±1; thus j) = 
(u(i),U(j)) in the fixed Cartesian reference frame. 

We recall that any product lattice, such as the D2Q9, is characterized by natural moments. 
For D2Q9, these natural moments are pMpq, where p = {f(ip)) is the density, and 

pMpq = P,9 6 {0,1,2}. (6) 

In the sequel we use the following linear combinations to represent natural moments 

Mqq^Ux = Mio, Uy = Mqi, T = M 20 + Mo2, N = M 20 - Mo 2, ^xy = Mu, 

Qxyy = Mi 2, Qyxx = M21, A = M22. 

These are interpreted as the normalization to the density (Mqo = 1), the flow velocity components 
{ux, Uy), the trace of the pressure tensor at unit density (T), the normal stress difference at unit 
density (N), and the off-diagonal component of the pressure tensor at unit density (H^jj,). The 
(linearly independent) third-order moments {Qxyy, Qyxx) and the fourth-order moment (A) lack 
a direct physical interpretation for incompressible flows. 

The b linearly independent moments serve as a different basis for the kinetic equations. It is 
clear that the choice of basis vectors is not unique. Another basis is given by the central moments 
of the form 


pMpq - Ux)^{vq^ Uy)'^f(iJ)), ( 8 ) 

where we have the following relation between the natural and the central moments 

Hxy “ Hajy UxUy, 

N = N+{ul-ul), 

T=f + u‘^, 

- - , - , - 2 (9) 

Qxyy ~ Qxyy ‘^UyA-xy ~ ^Ux^^ UxUy, 

Qyxx ~ Qyxx ‘^UxAxy + 2^y^ ^y^x, 

A — A 2 ^UxQxyy ^UxUyYl^y 2^ T “ 2^^^ ^y)^ ^x^y' 

We remark in passing that the mapping of natural moments onto central moments is nonlinear 
(it explicitly depends on the powers of the velocity components). 

It is important to note that the macroscopic equations are recovered by a projection of the 
kinetic system onto the lower order moments p,Ux,Uy. The higher order moments are thus in 
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our hands, in principle. However, they play an important role for the numerical stability of the 
scheme. Precisely this observation has been widely used to construct models with the goal to 
stabilize the LB scheme, among which are the MRT and KBC models. 

With the set of natural moments Q , populations are uniquely represented as follows (a, A = 

{- 1 , 1 }): 


/(o,o) = p{^-T + A), 

/(<T.O) = + ^) + - aQa:yy - A) , 

Ao.a) = 1^) + - ^Qyxx - a) , 

f ((7^X) ~ '^P A (^) (A)na,y + (jQxyy ^Qyxx ) • 


A similar representation can be written using the central moments by substituting eq. © in @. 

We group the population’s natural and central moment representations into the following 
functions for convenience: kt, ti{T), ni{N), Pi{I\.xy), qi{Qxyy,Qyxx), ai{A) and ki, ii{f), ni{N), 
Pi{Axy), qi{Qxyy,Qyxx), ai{A), respectively, where each of these groups also depends on the 
density and velocity, p,u, of the flow. ki(p,u) and ki{p,u) represent the kinematic part only. 
Thus, the populations are rewritten as 


fi — ki + + PiiAxy) + Qii^Qxyy^ Qyxx^ ^2(^)5 (H) 

shown here in the natural moment representation. With this formulation, the mirror state ([^ 
can be redefined by introducing relaxation parameters jk 

/r" = K + + (1 - lo)U{T)] + + (1 - jMN)] 

+ \l2Pi{Aby) + (1 - l2)Pi{Axy)\ + \l'iqi{Qbyy^Qbxx) + (1 “ ibqiiQxyyMyxx)] ( 12 ) 

+ [74ai(A®‘5) + (1 - 74)aj(A)]. 

We now need to find the optimal values for the relaxation parameters 7 ^. Let us first note that 
the relaxation parameters for the parts depending on the entries of the stress tensor, n^, Pi, must 
be set to 7 i = 72 = 2 in order to reproduce the Navier-Stokes equations, at least to second order. 
The relaxation for remaining moments can be chosen without affecting the hydrodynamic limit. 

A seemingly natural choice would be to set these relaxation parameters to 7 ^ = 1//3 (that is, 
7 i R, 1 for /? ->• 1) for i e {0,3,4} which implies that all higher order moments are brought to 
their equilibrium. This is essentially the idea of regularized lattice Boltzmann model uniiaiin]- 
However, in many benchmark simulations this choice does not lead to significantly better results 
as compared to the LBGK. MRT models on the other hand suggest highly optimized but fixed 
values for the relaxation of higher order moments. 

The KBC models take a different route by making the relaxation of higher order moments 
adapt to the flow and by letting local entropy decide about the corresponding relaxation values. 
We will discuss here KBC models with only one free relaxation parameter where the populations 
are represented as sum of three moment functions 


fi = h + Si + hi 


(13) 


where b (= kinematic part) depends only on the locally conserved fields, Si (= shear part) 
depends on the stress tensor 11 = Ei=i® "fi/i, and b (= higher-order moments) is a linear 
combination of the remaining higher-order moments. In the presentation (131 we essentially lump 
together all the higher order moments. Further extensions can be envisaged by splitting the 
higher order moments, hi into individual components in order to further improve the stability. 
However, we show in this paper that even with the suggested lumping of moments hi extremely 
stable LB models can be readily created which outperform MRT models. 
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Model s{-;p,u) 

h(--p,u) 

LBGK , A^, Tf Qxyy , Qyxx , ^ 

- 

KBC A n^y, N,f 

Qxyy-! Qyxx-! ^ 

KBC B n,,a,A 

^5 Qxyy! Qyxx-, 

KBC C n^y, N,T 

Qxyy! Qyxx! A 

KBC D n,,y,A 

r, Qxyy! Qyxx! 


Table I: Moment grouping for discussed models. 


With the representation (131, the mirror state is in a one-parameter form, 


/r'- = fc. + [2sr-s.] + [7/ir+(i-7)/^.], 


(14) 


where 7 is a relaxation parameter which is not yet specified. For 7 = 2 , KBC ( |14[ ) coincides with 
LBGK. For any 7 , the resulting LB model still recovers hydrodynamics with the same kinematic 
viscosity v (|§ which is demonstrated in section]^ 


III. MODEL DESCRIPTION 

Let us define the four KBC models A, B, C and D which we consider here. Models A and B are 
represented in the basis spanned by central moments while C and D are represented by natural 
moments. Models B and D, however, differ from A and C by including the moment function 
depending on the trace of the stress tensor, and t,(T), respectively, in the higher order part 

h. This eventually leads to different coefficients for the bulk viscosity. Table |I] summarizes the 
contribution for s and h for the different models. 


IV. ENTROPIC STABILIZER 


We will review the definition of the entropic stabilizer 7 given by m in the following. Let 
S{^) be the entropy of the post-collision state appearing on the right hand side of (Q, with the 
mirror state (141. We require the stabilizer 7 to correspond to the maximum of this function. 
Introducing deviations Asi = Si - and Ahi = hi - the condition for the critical point reads: 


^ Ahilnl 1 + 


(1 - Pj)Ahi - (2/3 - 1) As*' 

pq 

J i t 


= 0 . 


(15) 


Equation (151 suggests that among all non-equilibrium states with the fixed mirror values of the 


stress, = 2 s®‘^ - Si, we pick the one which maximizes the entropy. Note that 7 self-adapts 
to a value given by the maximum entropy condition at each grid point (15) and thus eliminates 


the need for tuning. We observe in all our simulation that the system entropy is monotonically 
growing over time. It is therefore conjectured that the second law of thermodynamics is fulfilled in 
practice and moreover, by providing a Lyapunov function, the entropy based relaxation contributes 
crucially to the stability of the scheme. 

An estimate for 7 in a closed form can be obtained by introducing the entropic scalar product 
{X\Y) in the 6 -dimensional vector space, 

b Y V. 

{X\Y) = Z^^ 

1=1 Ji 


( 16 ) 
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and expanding (151 to the first non-vanishing order in Asij and Ahij which yields 


* 1 /g l\ (As|Afe) 

^ P \ {Ah\Ah) 


(17) 


This estimate has proven to be sufficient for all practical purposes by stabilizing the scheme and 
monotonically incrementing the system entropy S. It is used for all the simulations presented in 
this work. 

Unlike any MRT, the relaxation parameter 7 is neither fixed a priori nor is it constant 
in space and time for KBC models. Therefore, it is instructive to obtain an estimate on 
the asymptotics for the different models. Let us expand (As|Aft,) in powers of velocity u, 
(As|A/i) = (As|A/i)*^°^ + (As|A/i)*^^^ + •••. At zeroth order this leads to 


(As^|A/r^)(°) = -ip(9A-l)(3r-2) 

(18) 

(Asb = 0 

(19) 

(Asc|A/rc)(°^ = -ip(9i- l)(3r- 2) 

(20) 

(Asb A/ib)^°^ = 0 

(21) 


where we have chosen to replace natural moments by central moments according to ([^. Note 
that the two models, B and D, which include the trace of the stress tensor in the higher-order 
part h do not contribute to the leading order. The next order contributions are 


(As^|A/i^)(i) = 0 (22) 

(Asb|A/ib)^^^ = 0 (23) 

{AsclAhc)^^'’ = lp{Qxyy{{2 + 3iV - 3f)ux - l2p^yUy) - Qyxx{^‘2PxyUx + (-2 + 3N + 3f)uy)) 

(24) 

{Asj^\Ahj)'j^ ^ ^ p(^~NQxyy'^x ^PxyQyxx^x '^PxyQxyy^y P^Qyxx^y^ • (^^) 

We see that KBC B is the only model among the four that does not contribute to either zeroth 
or first order. This leads to the conjecture that, when /3 -> 1, 7 will be close to 1 for the KBC B 
model, unlike the other KBC models. This is clearly confirmed by our simulation results (see, e.g. 
figs. 0 and . The average value of the entropic stabilizer 7 fs 1 is found in our simulations for a 
range of Reynolds numbers and resolutions. Note that by fixing 7 = 1//3 the KBC models coincide 
with the regularized LB model (more specifically, with a realization of the regularized LB in the 
according moment basis). Thus, the above theoretical derivation of the entropic stabilizer for 
KBC model B produced the empirical regularized LB in the central moment basis. 


V. HYDRODYNAMIC LIMIT OF KBC MODELS 


Let us derive the hydrodynamic limit of the general kin etic equation with KBC-type mirror 
state (141. We start by rewriting eq. (Q and eq. ((T4| as 


/' = /, +2/3 (/f^=-/,). 

with the generalized equilibrium [3111121] of the form. 


(26) 


_^GE 


/r+i( 7 - 2 )(/ir-^o- 


(27) 
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Model (siVixVix) (siViyViy) 

i^SiVixViy 

) {hiVix^ix 

) {hiViyViy 

) {hiVixViy} 

KBC A \p{f + N) \p{T-N) 

pTVxy 

0 

0 

0 

KBC B \pN -\pN 

pWxy 

ipT 

\pT 

0 

KBC C \p{T + N) \p{T-N) 

P^xy 

0 

0 

0 

KBC D \pN -\pN 

pYVxy 

hpT 

\pT 

0 


Table II: Second order moments for functions s and h in KBC. 


In the following derivation, Einstein’s summation convention is applied for all subscript indices 
except for i where the explicit notation (...) is used. It is useful to compute the second and third 
order equilibrium moment until second order in velocity beforehand 


= (fi'^V.aVip) = pclSa)3 + (28) 


Due to local conservation laws and as a direct consequence of the construction of the moment 
groups k, s and h (see eq. (131 and table we can immediately state the following relations for 
the zeroth and first order moments 


\_Pj P^a\ 1 


(30) 

(31) 


While these relations hold for all four KBC models discussed here, they depart from each other 
in the higher order moments. Table ^ shows the second order moments for moment functions s 
and h, respectively. Note that all the higher order moments for the kinematic part k vanish. 

After the previous preliminary considerations let us expand the left hand side of equation (26 1 
using a Taylor series to second order 

^St{dt + daVia) + ^{dt + daVia){dt + ^/S^^i/?)] fi = 2/3 [/°‘* - fi+ |( 7 - 2)(ft,®‘* - hi)] . (32) 


By introducing a characteristic time scale of the flow, 0, we can rewrite ( |32[ ) in a non-dimensional 
form using reduced variables t' = tjQ, v[ = Vijc and x' = xl{cQ), where c = 1. After introduction 
of the parameter e = St/Q and dropping the primes to simplify notation we get 

[e(i9t + daVia) + Yidt + daVia){dt + dfjVip)'^ fi = 2j3 [/®‘^ - /^ + §( 7 - 2)(/i-‘’ - hi)]. (33) 

By exploiting the smallness of e we can perform a multiscale expansion of the time derivative 
operator, the populations and their decomposition into s and h parts until second order. 



(34) 


(35) 

Si = + ■" 

(36) 

h, = h^KePp^e^hpK.... 

(37) 


Inserting eqs. (34) to (37) into eq. (33) we can analyze the terms corresponding to orders e°, 
and e^. The zeroth order terms lead to 


0 = 2/3[/r-/f' + |(7-2)(/ir-/ir^)], 


( 38 ) 










which implies 




Local conservation laws dictate the relations {fi{^,Via}) = and 

{hi{l,Via}) = {h‘l‘^{l,Via}) = 0 which yield the solvability conditions 





= 0 , 
= 0 , 


^ = -2f3 [/f) + i(7 - 2)h^^], 


( 39 ) 


(40) 

(41) 


(42) 


from which we can recover the hydrodynamic equations of mass and momentum to first order by 
taking the zeroth and first order moment of (421 and using condititions (401 and (411 and the 
definition of 11®'!, 

ap 


(43) 

(44) 


p = -da{pUa) 

= ^^Uadp{pup) - 

Collecting terms of order results in 

[af ^ + daVia){dl^^ + dpvip)^ Z®'^ + = -2/3 [/f ^ + §(7 - 2)hf . (45) 


From eq. (421 we get the relation 

+ 1(2 - 


which, once inserted in eq. ( |45| ), leads to 

[>9^^ + (i - + daVio,){dl^^ + dpvip)^ Z®'’ + + 9«?;*a)(2 - 7)hf ^ 

= -2^[zf) + |(7-2)/if)]. 


(46) 


(47) 


Taking the zeroth order moment thereof and making use of eqs. ( |40[ ), (41), ( |43| and ( |44[ ) yields 
the vanishing second order contribution to the continuity equation 


afV-o. 


(48) 


The first order moment, however, will render different outcomes for the four KBC models 
depending on the second order moment of h, 

= ^ (^ - I) 9/3 + ^dp [(7- 2 )! (/if ^u,„^;i/3)]. ( 49 ) 

For models A and C the last term on the right hand side vanishes according to table [Tl] and we get 

= 19^ [( A_ _ 1 ) (af )n®q + , (50) 
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while for models B and D one must first analyze the second moment of h. As only the trace 
ViaVia^ will contribute we can get the following relation using eq. (46 1 , condition = 


and table [fl] 


ViaVip'^ = I Vi^Vi^^ 5ap = “ (^ (sf ^ + d^Via) SajS 


(51) 


= -—+8 « 

2'y(3\^t ^oc(3’ 

When inserted in eq. ( |48[ ), this gives the second order contribution to the momentum equation 

For completeness let us now compute the first order terms of the pressure tensor for models A, 
C and B, D which yield, by substituting the equilibrium values eqs. (28) and (29) and using the 
first order results eqs. (43) and (44) 


(52) 


nf .^2 = 




(53) 


and 


nia = [daUp + dpUa - ^ [^d^u^Sa/s] , 


(54) 


respectively. 

By summing up the contributions from first and second order in e, using the eqs. (28 1 , (291, (431, 
(44) and (48) and reintroducing dimensional variables we recover the isothermal Navier-Stokes 
equations at reference temperature Tq = 


dtP = -daipua), (55) 

with the kinematic (shear) and bulk viscosity coefficients 

2 / ^ ^ \ \ v models A and C 

^ = 2 ). ij models B and D. 

Note that for the quasi-incompressible LB method the bulk viscosity term is small and can be 
considered as an artefact of the numerical method. We also point out that the bulk viscosity ^ 
for KBC B and D depends on time and space along with the stabilizer 7 . 


VI. GREEN-TAYLOR VORTEX FLOW 

For all of the flows considered in this paper periodic boundary conditions are applied 
in order to separate the accuracy of the scheme from influences of boundary conditions. 
The KBC scheme is validated for the Green-Taylor vortex flow in a first numerical ex¬ 
ample. The analytical solution for the Green-Taylor vortex flow is given by u{x, y, t) = 
V X [(mo/A: 2 ) cos(fcia;) cos(fc 2 y) exp(-:/(fcj + fc|)t)] where we have chosen A:i = 1, fc 2 = 4 and the 
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Figure 1: Convergence rate for Green-Taylor vortex at < = tc- (a) Re = 100, uq = 0.03, (b) 
Re = 1000,uo = 0.04. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), KBC D 
(□), ELBM (o), second order convergence (fine dotted). 


pressure po = P<?s is initialized using po = 1- The populations were initialized using Grad’s 
approximation m 


f* {P,U,n) = Wi [p+ + ^(Ilap - pclSap)iViaVii3 - C^l^a/?)] 


(58) 


where the pressure tensor was taken in the form (see eqs. (28) and 

is confined in 0 < a:,?/ < 27r covered by a mesh of fV x iV lattice nofes, the 
is defined as Re = uqN jv and the decay half-time of the flow is given by tc = 
lattice time steps. 

The Reynolds number was set to Re = 100 in a first experiment and Re = 
simulation while the initial velocity was mq = 0.03 and uq = 0.04, respectively, 
varied between {64,128,256} in the former and between {62,128,256,512} 
KBC models A-D, LBGK and ELBM were run and compared to the analytic 
shows the convergence for the relative error at time t = tc- Second order rate 
for all models, moreover, all the models are performing almost identically. 


]53)). The domain 
Reynolds number 
= ln2/[i/(fc2 + fc|)] 

1000 in a second 
Resolution N was 
in the latter case, 
solution. Figure [T] 
is clearly observed 


VII. DOUBLY PERIODIC SHEAR LAYER 


To probe the KBC models for their performance in under-resolved simulations of smooth flows 
with sharp features the doubly periodic double shear layer with a slight perturbation studied 
extensively in m was used as a benchmark. Initial conditions are given by 


Uy 


f Mo tanh (k (p/V - 0.25)), y < fV/2, 
\ Mq tanh (k (0.75 - yIN)) , y > N12, 

Suo sin (27r (x/N + 0.25)). 


Here N is the number of grid points in both x and y directions while periodic boundary conditions 
are applied in both directions. Grad’s approximation (581 was used to initialize the flow field 
while initial density was set to unity. The parameter k controls the width of the shear layer 
while (5 is a small perturbation of the velocity in y-direction which initiates a Kelvin-Helmholtz 
instability causing the roll up of the anti-parallel shear layers. Uq is the initial magnitude of 
the a;-velocity while the Reynolds number is defined as Re = uoNjiy and the turnover time is 
tc = NIuq lattice time steps. 










11 



Figure 2: Vorticity field at t = tc for k = 80, uq = 0.04. Columns: LBGK, KBC A, KBC B, KBC 
C, KBC D, ELBM, respectively, rows: Resolution N = 128,256,512. 


In [28] it is demonstrated that all the numerical methods investigated therein produce spurious 
additional vortex roll-ups as a consequence of under-resolution. Effectively no convergence could 
be reported until the resolution was sufficiently high for the additional vortices to disappear. 

We pose the question whether the different relaxation for the non-hydrodynamic higher order 
moments are advantageous for the performance in under-resolved cases. To this end let us 
consider a thin shear layer case with k = 80, Re = SO'OOO and uq = 0.04. We compare KBC models 
A-D, LBGK and ELBM. Figure shows the vorticity field at t = for the six schemes under 
consideration for different resolutions N = {128,256,512}. LBGK becomes unstable even before 
t = tc is reached (see also fig. [^a) and b)) for N = 128. For N = 128, model C and ELBM clearly 
show formation of additional roll-ups whereas model D produces comparatively small instabilities. 
Models A and B (central moments) capture the flow features quite accurately while model B seems 
to perform slightly better. For the next higher resolution under consideration, N = 256, LBGK 
survives but features two small additional vortices while the other models capture the main flow 
features well. For the largest resolution, N = 512, the models are essentially indistinguishable. 

In summary, KBC B performs qualitatively better in the under-resolved situation than the 
other models while for the still slightly under-resolved case, N = 256, the KBC models and ELBM 
give comparably good results while LBGK still features spurious vortices at this resolution. 

These findings are also reflected in figure ^ where the second order convergence is reached for 
all models after N = 256. KBC models B and A clearly outperform the other schemes in the 
under-resolved cases. 

Let us now consider the energy and enstrophy decay, eq. where we report both the mean 
and the fluctuations (RMS) over time. We first remark that the methods converge to each 
other for N = 256 while there is also evidence that the simulation is resolved as the statistics do 
not change for the next higher resolution, N = 512. The energy decay is rather similar for all 
the models across all resolutions, indicated by both mean and standard deviation, except for 
ELBM which shows slightly different results for the mean and fluctuations at N < 128. The mean 
enstrophy and fluctuations over time are clearly better captured for the KBC models A, B and D 
in the under-resolved situation compared to ELBM and KBC C. This is also in accordance with 
the visual impression of the vorticity structure (fig. |2 I. In summary, the low order statistics seem 
not affected by the KBC treatment of the higher order moments. 




12 




t/tc 



t!tn 



t/tc 


Figure 3: Evolution of kinetic energy and enstrophy (mean: left axis and large symbols, standard 
deviation: right axis and smaller symbols). Resolution: a) and b) iV = 128, c) and d) iV = 256, 
respectively. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), KBC D (□), ELBM 
(o). Results represented by symbols have been subsampled for clarity. 



Figure 4: Convergence rate for doubly periodic shear layer at t = and 
Re = 30000,K = 80,Mo = 0.04. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), 
KBC D (□), ELBM (o), second order convergence (fine dotted). Error with respect to reference 

solution (LBGK at = 2048 resolution). 
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Figure 5: Vorticity field for decaying two-dimensional turbulence Re = 13'134, N = 1024, for 
LBGK (first row) and KBC B (second row) and times tjte = {10,20,40,60,80,100} from left to 

right. 


VIII. DECAYING TWO-DIMENSIONAL TURBULENCE 

The third and final numerical example considered in this paper is the simulation of a turbulent, 
albeit two-dimensional, flow starting from a random initial condition and decaying with time. 

Decaying two-dimensional turbulence is characterized by the formation of vortices in the 
early stage (vortex generation period) which leads to spatially separated coherent structures 
(which account for the vorticity extrema) which typically have long lifetimes compared to the 
eddy turnover time and undergo passive advection and vortex-vortex interaction |26| . These 
vortices can persist, grow over time as they merge with weaker structures of same-sign vorticity 
and influence the whole field [S] . A large amount of the enstrophy is concentrated within the 
large-scale vortices that decay slower than the background vorticity field between the vortices 
m- The total energy is roughly constant while enstrophy is decaying. According to m the 
enstrophy follows a direct cascade from large to small scales, much alike the energy cascade in 
three-dimensional turbulence, while the energy shows an inverse cascade from small to large 
scales. Classical Kolmogorov-Batchelor scaling theory predicts a slope of k~^ and k~^ for the 
energy (E) and enstrophy (Z) spectrum, respectively, where k is the wave-vector magnitude. 

In general, for the simulation of fluid turbulence it useful to study the following questions: 

1. is the dynamics of a fully developed turbulent flow accurately captured (i.e. initial vortex 
formation, emergence of coherent structures, vortex-vortex interactions, decay of vortex 
density)? 

2. are near grid-scale structures with large gradients (small vortices) well represented for 
sufficiently high Reynolds numbers? 

3. is the stability affected by the these lattice-scale structures? 

4. are low-order statistics well represented and is the numerical scheme correctly modelling 
the physical dissipation (i.e. decay of enstrophy, scaling laws)? 

5. how good is the performance for very large Reynolds numbers in an under-resolved simula¬ 
tion? 

The initial conditions for all subsequent simulations are given by constructing a zero-mean 
Gaussian random field in Fourier-space with random Fourier-phases and amplitudes proportional 
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Re= 13'134 


Re= 1.5-10® 

Re = 1.6-10® 


N 256 

512 

1024 

1024 2048 4096 

1024 2048 

4096 

{B(0)> 1.645 ■10"“ 
(Z{Q)) 8.455-10"® 
V 3.533 ■ 10"“ 
te 344 

1.645 ■ 10"“ 
2.137-10"® 
7.066 ■ 10"“ 
684 

1.645 ■ 10"“ 
5.358 ■ 10"^ 
1.413 ■ 10"® 
1366 

1.209-10"“ 1.209-10"“ 1.209-10"“ 
3.815-10"® 9.602-10"® 2.404-10"® 
1.062-10"“ 2.123-10"“ 4.246-10"“ 
512 1021 2039 

1.209-10"“ 1.209-10"“ 
3.815-10"® 9.602-10"® 
9.952-10"® 1.990-10"® 
512 1021 

1.209 ■ 10"“ 
2.404 ■ 10"® 
3.981 ■ 10"® 
2039 


Table III: Characteristics for two-dimensional turbulence simulations (lattice units). 


to the prescribed spectral density of the stream function 4'(A:) = k~^Z{k) = k~^E{k) from which 
an incompressible velocity field for a domain of N x N lattice nodes is obtained. The energy 
spectral density function is given by 

E{k) = Cok^[l + {k/ko)^^^Y^ (59) 

where Cq is a normalization constant and the parameters A = 6 and B = 17 such that the energy 
spectrum is narrow banded and reasonably peaked at small wave numbers [3]. 

We consider three groups of numerical experiments. The first simulations are carried out at a 
moderate Reynolds number Re = 13'134 where the energy spectrum was peaked at fcg = 9, the 
second group was run at Re = 1.5 • 10® with a different random initial velocity field, ko = 30, and 
the last group was simulated at a very high Reynolds number Re = 1.6 • 10® with the same initial 






Figure 6: Evolution of kinetic energy and enstrophy for Re = 13'134. Mean: left axis and large 
symbols, standard deviation: right axis and smaller symbols. Resolution N = 256,1024 from top 
to bottom, respectively. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), KBC D 
(□), ELBM (o). Results represented by symbols have been subsampled for clarity. 
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Figure 7: Evolution of mean stabilizer 7 , mean entropy (right column, left y-axis, large symbols 
and right y-axis, small symbols, respectively) and palinstrophy (left column, mean: left y-axis, 
large symbols, standard deviation: right y-axis, small symbols) for Re = 13'134. Resolution 
N = 256,1024 from top to bottom, respectively. LBGK (solid), KBC A (dashed), KBC B (x), 
KBC C (dotted), KBC D (□), ELBM (o). Results represented by symbols have been subsampled 

for clarity. 




Figure 8 : Energy (a) and enstrophy spectra (b) for two-dimensional turbulence at Re = 13'134 for 
N = 1024 at time tjU = 50. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), KBC 
D (□), ELBM (o). Results represented by symbols have been subsampled for clarity. 
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Figure 9: Evolution of kinetic energy (a), enstrophy (b), palinstrophy (c) and stabilizer 7 (d) for 
Re = 1.5 • 10® at resolution N = 4096. Mean: left y-axis and large symbols, standard deviation: 
right y-axis and smaller symbols. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), 
KBC D (□), ELBM (o). Results represented by symbols have been subsampled for clarity. 


field. Reynolds number is defined as Re = N\J{2E)li>, where E is the mean initial kinetic energy. 
A rough estimate of the eddy turnover time is given by te ~ [3]. All the simulations were 

run for 100 tg in order to observe both vortex formation, merging and decay. The characteristic 
figures for the initial conditions are summarized in table [ml As in all our simulations. Grad’s 
approximation (58 1 was used to initialize the populations, while here the gradients of velocity 
were estimated by central differences from the given initial field. Note that LBGK could not cope 
with the under-resolved cases N = 1024,2048 at Re = 1.6 • 10® and “crashed” due to numerical 
instabilities where the other five methods run trouble-free. 


Let us first consider the low Reynolds number case. Figure shows a comparison of the 
vorticity field for LBGK and KBG B and N = 1024 at different times. The first column shows the 
vorticity structures at the point of maximum turbulence activity indicated by the palinstrophy 
evolution (see fig. [^. The vortices have been formed and coherent structures appear in the next 
shown time instance which interact with each other. The number of vortices is clearly decaying 
when comparing the last column of fig. |^with the earlier time instances. It is striking that all 
the models show almost the same dynamics (for brevity fig. shows only LBGK and KBC B). 
Except for the last time point the plots are visually hardly discriminable. Even after a long time, 
tjte = 100 , the vorticity structures are still comparable. 

The decay of enstrophy was measured and figs, [^b) and d) show the expected exponential 
decay for resolutions N = 256 and N = 1024. It is apparent that the evolution of mean enstrophy 
is the same for all models and is almost identical among the two resolutions. The fluctuations are 
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Figure 10: Energy (a) and enstrophy spectra (b) for two-dimensional turbulence at Re = 1.5 • 10® 
for N = 4096 at time tjte = 50. LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), 
KBC D (□), ELBM (o). Results represented by symbols have been subsampled for clarity. 


largely the same, although one can observe a slight flattening at later times for the KBC models 
compared to ELBM and LBGK for N = 256. The evolution of energy, figs. §a) and c), shows 
a similar tendency; the models coincide in the mean but differ in the fluctuations for the lower 
resolution. As all the models are close to each other at N = 1024 and the means of enstrophy and 
energy are not changing compared to N = 512 we consider this highest resolution as resolved. 

Evidence for this classification is also gathered from figure a) and c), where the mean 
palinstrophy and its fluctuations match for all models at the highest resolution. As stated 
earlier, at tjtf, « 10 we observe a peak in palinstrophy which indicates a state of high turbulence 
intensity. Note that the maximum value is slightly better captured by ELBM and LBGK in the 
low resolution case. 

For all the measured low-order statistical moments we have seen almost identical values, at 
least in the resolved case, and largely identical mean statistics overall. It is interesting to see, 
however, that the KBC models differ quite significantly among each other with respect to the 
evolution of the stabilizer 7 , see fig. [^b) an d). Especially, KBC B is fundamentally different 
from the other KBC models by 7 staying close to 1 all the time. This is consistent with the 
“quasi-orthogonality” of the A/i and As decomposition of this model demonstrated above (see 
eqs. (191 and ([^). Nevertheless, the overall production of entropy is nearly identical for all the 
considered models. 

The scaling of spectral energy and enstrophy density is shown in fig. for N = 1024 and 
tjte = 50. Due to the moderate Reynolds number the slopes are not expected to match with 
the theoretical prediction for high Reynolds numbers, however, we are not able to distinguish 
between the models. Thus we can conclude that for the moderate Reynolds number all methods 
give almost identical results and show the same dynamics. 

The next simulation was carried out with a higher Reynolds number. Re = 1.5 • 10®, but on 
a larger grid. The number of initial vortices is much higher due to a different random initial 
condition. Here we report the results for the highest resolution, N = 4096. For the lower 
resolutions, the results are similar, albeit with slightly more variance among the models with 
respect to the fluctuations of energy. 


According to figs. [^and[^ where evolution of kinetic energy, enstrophy and palinstrophy as 
well as spectral density of energy and enstrophy are reported, we observe that the six models 
behave almost identically. Note that for this Reynolds number LBGK did not encounter any 
numerical instabilities. Due to the similarities among the models and the fact that the results for 
N = 2048 are not significantly different from the largest resolution, we conclude that for N = 4096 
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Figure 11: Vorticity field for decaying two-dimensional turbulence, Re = 1.6-10®, for different 
ELBM (a) and KBC B (b) for time tjte = 50. 


the flow is essentially resolved. Note the differences among the models in figure [^d) for the 
stabilizer 7 which is in accordance with the results from the lower Reynolds number. Although, 
the Reynolds number here is one order of magnitude higher than before, the scaling of energy 
and enstrophy is still too steep compared to the theoretical slope. 

In order to verify that the discussed models can achieve the proper scaling laws for very high 
Reynolds numbers, we conducted a simulation at Re = 1.6 • 10®. Let us first remark that for 
this highly turbulent regime, LBGK was not able to run with N = 1024 and N = 2048. There 
is evidence that all of the considered grids do not fully resolve the flow as the mean statistical 
quantities are still slightly different among the two highest resolutions, N = 2048,4096, so that the 
six models are affected differently by the lack of resolution. It is thus interesting to see whether 
the dissipation is affected by the KBC models in in the presence of under-resolution. 

For this matter let us compare the ELBM model and KBC B for the time t/fg = 50 and 
N = 4096. Figure [m shows the vorticity field, accordingly. Note that the all the models produce 
approximately the same vortex structures up to t/te ~ 20. Although, one can still see similarities 
of the structures at tjte = 50, the two models produce distinctly different pictures. The number of 
vortices, roughly estimated by the number of vorticity patches exceeding two times the standard 
deviation of vorticity, is clearly different: ELBM accounts for 3440 whereas KBC B shows 1587 
vortices. While visually the number of larger vortices seems comparable, this difference must 
stem from the very small structures. 

These findings are also consistent with the energy and enstrophy spectra depicted in fig 
While ELBM keeps more energy and enstrophy in the large wave numbers, KBC models smoothly 
fall off. At lower resolution, ELBM shows a bump near the largest wave numbers (see fig. [T^ 
which was observed in other simulations as well. This is conjectured to be the effect of the built-in 
subgid model established through the fluctuating effective viscosity. It can be observed, however, 
that all the models capture the theoretical slope in a range of wave numbers. 

On the other hand, the enstrophy evolution, depicted in fig. 12 b), shows almost identical 
average dissipation for all the models at A^ = 4096. This indicates that despite the difference in 
the spectra KBC models do not introduce a significantly higher dissipation, however, it seems 
that the flux of energy to larger scales is more dominant than in the case of ELBM. Fig. 12 a) 
shows the enstrophy decay for the resolution N = 1024. Here, the ELBM method clearly shows a 
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Figure 12: Evolution of enstrophy for Re = 1.6 • 10® at resolution N = 1024 (a) and N = 4096 (b). 

Mean: left y-axis and large symbols, standard deviation: right y-axis and smaller symbols. 
LBGK (solid), KBC A (dashed), KBC B (x), KBC C (dotted), KBC D (□), ELBM (o). Results 
represented by symbols have been subsampled for clarity. In (a) the dark dashed line represents 
the corresponding values at = 4096 for LBGK. 



Figure 13: Energy (a) - (c) and enstrophy spectra (d) - (f),for two-dimensional turbulence at 
Re = 1.6 • 10® for N = 1024,2048,4096 (top to bottom) at times t/tg = 50. LBGK (solid), KBG A 
(dashed), KBC B (x), KBC C (dotted), KBC D (□), ELBM (o). Results represented by 
symbols have been subsampled for clarity. 


slower decay than KBC. When considering the corresponding curve for N = 4096 and LBGK as a 
reference (dark dashed line) it is apparent that none of the models capture the expected rate, 
however, ELBM slightly under predicts the decay at later times while the KBC models show in 
general more dissipation (see also fig. [T^ d) and e)). In contrast to the smallest scales, the large 
and moderately small scales seem to be predicted well by all methods. 
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Stability and Performance 

In all the simulations reported in this paper as well as in preliminary tests and the simulations 
published in m we did not encounter a single case of numerical instability using the KBC models, 
even when running very under-resolved simulations for high Reynolds numbers. The ELBM 
method is non-linearly stable and thus it is not surprising that it runs trouble-free for all the 
set-ups. LBGK on the other hand failed numerically in presence of under-resolution. 

The KBC models introduce additional computational overhead in order to compute additional 
moments and the estimate for the stabilizer 7 which accounts for not more than a factor of 2-2.5 
for both two and three dimensions. 


IX. CONCLUSIONS 


In this work we studied four variations of entropy based multi relaxation models of the recently 
introduced KBC family. We reviewed the details of the entropic stabilization and described the 
four models in detail. The recovery of the Navier-Stokes equations was demonstrated to second 
order for all the KBC models. A detailed comparison with LBCK and ELBM was carried out at 
various grid resolutions and Reynolds numbers for different two-dimensional flows. Second order 
rate of convergence is numerically confirmed for all the models studied herein. 

It must be stressed that the entropic models, KBC and ELBM, were stable (in contrast to 
LBCK) for all the considered cases here; despite of under-resolution and high Reynolds numbers 
(e.g. Re = 1.6 • 10® for a grid of 1024 x 1024). 

The accuracy of various models was discussed, in particular, for under-resolved simulations. 
We first remark, that all the models converge to the LBCK solution for resolved flows and that, 
in general, the difference among the models under consideration are small. At moderately high 
Reynolds number, the KBC models perform better than LBCK, which suffers from numerical 
instabilities, and ELBM which produces spurious vortices in the case of the double shear layer. 

The simulation of two-dimensional turbulence was chosen as a benchmark to assess various 
statistics for high Reynolds numbers. While covering a range of different resolutions and Reynolds 
numbers, the KBC models were shown to capture the expected scaling laws for energy and 
enstrophy spectra as well as ELBM (and LBCK for sufficiently large resolutions). There is 
indication that the KBC models produce less small structures than ELBM (and LBCK) for 
the under-resolved cases, however, the decay of enstrophy is only slightly accelerated in the 
under-resolved case for high Reynolds numbers. On the other hand ELBM tends to somewhat 
amplify the appearance of small structures which lead to a slight over-representation of enstrophy 
(and energy) content at large wave numbers but only for very coarse resolutions. In summary, all 
the KBC models considered here have the correct limit of LBCK for resolved simulations, but 
more importantly, they capture all low order statistics such as averages and fluctuations of kinetic 
energy, enstrophy and palinstrophy as well as the spectral densities for energy and enstrophy 
extremely well despite severe under-resolution. 

Among the KBC variations, model B was demonstrated to be more accurate in the under¬ 
resolved shear layer case. The significant difference with respect to the evolution of the stabilizer 
7 for model B can be explained by the orthogonality of the entropic scalar product (As|Ah) in 
the leading and first order terms in velocity powers. 

In general, we show that by keeping the kinematic (shear) viscosity coefficient constant (in 
contrast to ELBM) the presented method is extremely stable and produces accurate results in 
presence of under-resolution (similar to ELBM). Minor differences in performance among the 
different KBC versions are observed for different simulations, however, all KBC models are much 
more stable than LBCK. 
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It has been demonstrated in [22] that also for three dimensional flows in presence of complex 
walls low-order statistics can be captured well using KBC. In a further publication we will address 
the issue of boundary conditions for KBC models in both two and three dimensions. 
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